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Abstract We review and introduce a generalized reaction-diffusion approach to epi- 
demic spreading in a metapopulation modeled as a complex network. The metapopula- 
tion consists of susceptible and infected individuals that are grouped in subpopulations 
symbolising cities and villages that are coupled by human travel in a transportation 
network. By analytic methods and numerical simulations we calculate the fraction of 
infected people in the metaopoluation in the long time limit, as well as the relevant 
parameters characterising the epidemic threshold that separates an epidemic from a 
non-epidemic phase. Within this model, we investigate the effect of a heterogeneous 
network topology and a heterogeneous subpopulation size distribution. Such a system 
is suited for epidemic modeling where small villages and big cities exist simulta- 
neously in the metapopulation. We find that the heterogeneous conditions cause the 
epidemic threshold to be a non-trivial function of the reaction rates (local parame- 
ters), the network's topology (global parameters) and the cross-over population size 
that separates "village dynamics" from "city dynamics". 

Keywords Epidemic spreading • Diffusion • Complex networks • Heterogenous 
dynamics 



1 Introduction 

Since the beginning of man, infectious diseases have played an important role in how 
we interact, migrate, and form societies. Diseases that are highly infectious, can soon 
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develop into an epidemic. A "global epidemic" — often referred to as a pandemic — 
is a disease that cause infection in a significant portion of a population over a large 
geographical area. Fortunately, most of the common infectious diseases are not lethal, 
but the ones that are can result in huge numbers of deaths, often significantly reducing, 
and even eradicating, the population in certain areas. 

Historically, the world has experienced several pandemics, often with high death 
tolls; The plague of Justinian (541-542AD) struck the Eastern Roman Empire and is 
estimated to have had a mortality of 50% to 60%; the Black Death of the 14th century 
reduced the world population from an estimated 450 millions to between 350 and 375 
millions; the tuberculosis pandemic of the 19th century killed about 25% of the adult 
European population; and the Spanish flu (influenza) of 1918 is estimated to have 
caused the death of between 25 to 50 million people. More recently, the swine flu or 
the HlNl influenza (a new variant of the Spanish flu) killed about 300,000 people 
within the first 12 months in 2009 1 1 1 . In 2002, Hong Kong experienced an outbreak of 
the severe acute respiratory syndrome (SARS) which, according to the World Health 
Organization, nearly became a pandemic with 8,422 cases and 916 deaths worldwide 
between November 2002 and July 2003 (10.9% fatahty) |2|. SARS spread within 
weeks from Hong Kong to 37 other countries. 

Needless to say, there is a huge interest in epidemic outbreaks, spreading and con- 
tainment which has attracted scientists from a number of different disciplines. Early 
pioneers in infectious disease modeling were William Hamer and Ronald Ross (No- 
bel laureate in medicine 1902) who in the early twentieth century applied the law of 
mass action to explain epidemic behavior of common infectious fevers of childhood 
and malaria, respectively |3 ,4|. The early work by Kermack and MacKendrick (5l| on 
disease dissemination in a homogeneous population also deserves acknowledgment 
in this context. Since then numerous works have appeared of increasing complexity 
and detail. One set of models concern spatial dynamics dealing with epidemic spread- 
ing over geographic areas |6|. Examples range from forest-fire models 17) to spatial 
immunization patterns where multiple diseases compete for their hosts fSl, and ex- 
plicit simulation of 85 million people including households, schools and workplaces 
in Thailand Q, a region of high interest due to the avian H5N1 epidemic. 

Another class is metapopulation models. Such models are based on the obser- 
vation that large human populations are typically divided into a network of smaller 
subpopulations such as cities or villages. Crudely speaking, a metapopulation con- 
sists of a group of connected subpopulations. Between the subpopulations there is an 
ongoing exchange of people, some of which are carriers of a disease and may infect 
previously unaffected subpopulations p0|{r4[ . The dynamics in each "node" is often 
assumed to be governed by mass-action kinetics, or sometimes stochastic (relevant 
for very low population densities). Obviously, the way in which the subpopulations 
are interconnected has a significant effect on how the disease spreads through the 
metapopulation. Particular emphasis has therefore, on good grounds, been paid to the 
impact of airline travel using explicitly, or features of, the world's airline transporta- 
tion network p5[|T9 |. Recent efforts also try to improve the description of human 
mobility patterns such as recalling individuals' geographic origins |20||21) . 

A third model class, albeit related to the former, is network models that account 
for how humans actually meet p2] - p4| . In metapopulation models it is often assumed 
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that the population within a single subpopulation is well mixed, which means that 
every individual is equally likely to spread the disease to any other member of the 
population. This is of course not generally true. For example, in large cities any 
individual only meets a very small fraction of the total population. There is, however, 
frequent encounters between those belonging to the same social network. 

There are also interesting works focusing on the situation where humans are not 



the main carriers of the disease. Examples include Bubonic plague carried by rats |25 1 
and water-borne cholera epidemic in a South African river network p6| . 

In this paper, we use a metapopulation description to model disease spreading in 
a network of cites and villages (nodes) of varying population sizes which are inter- 
connected via a transportation network. The infection model we assume is the well- 
established susceptible-infected-susceptible model (Sec.[2]l where individuals can be 
in either a susceptible or an infected state and going from one state to the other is 
controlled by rates. For an isolated city or village, the fraction of infected individuals 
in the population in the final (stationary) state is fully determined by these rates. How- 
ever, when the cities and villages in addition are exchanging people with each other, 
for instance on a transport network, the situation may be much more involved. Under 
such conditions, the final state is not only determined by the rates, but also depends 
on structural features of the network 1 10, 1 1 1. 

The final state of the metapopulation can be either epidemic or non-epidemic, 
meaning that a fraction of the population is infected (or not). The separation between 
these two phases (in parameter space) is known as the epidemic threshold. In this study 
we go beyond previous works and address the importance of city-size heterogeneity 
on the epidemic threshold. The size of the city or village is an important aspect because 
it sets the limit on how large fraction of the population a single individual actually 
meet within a given time period. The well-mixed approximation works poorly for a 
large city like New York whereas for a small countryside village it could be adequate. 
We therefore let the infection dynamics at each node change depending on how many 
individuals are present at any given time. In previous models, the network's topology 
only influenced how people, and thus the disease, spread between the nodes. Here we 
in addition let the dynamics adapt to the local instantaneous population density, which 
means that the network structure influence the disease dynamics in two ways. First, 
how individuals move, and second, how they interact on the nodes. The latter part is in 
practise achieved by introducing a population dependent infection rate that interpolates 
between the two limiting cases of "the big city" for which the infection rate is dropping 
with increasing population, and "the small village" where the infection rate remains 
independent of population size. Within this setting, we derive new analytic expressions 
for the epidemic threshold of a general network which interpolates between previously 
known cases. 



This paper is organized as follows; Section[2]is devoted to detailing the infection 
model that will be assumed. Here also the behavior of this model in homogeneous 
space is discussed. We then present how we model the transportation of individuals 
between the cities and villages (Sec.[3]l before we in Sec.|4]present a unified reaction- 
diffusion model that combines the two phenomena from the previous two sections. 
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Section|4]also presents analytic results for the epidemic threshold and compare them 
to simulations. The conclusions from this study are presented in Sec. [5] 

2 The susceptible-infected-susceptible epidemic model 

The model that we will focus on in this work is the classic susceptible-infected- 
susceptible (SIS) epidemic model [27 J . This model is chosen since it is relevant for 
an important group of infectious diseases. Moreover, it has been studied extensively 
in the literature, and it can often be given an analytic treatment. 

In the standard SIS model there exists two types of individuals — the individu- 
als that are susceptible to an infection (labeled S below), and those individuals that 
already are infected and can infect others (labeled I). If an infected individual meets 
an uninfected individual, he will transmit the disease at rate /3. On the other hand, an 
infected individual can recover from the illness at rate fji and thereby again becoming 
susceptible to the disease. Formally, the reactions of the SIS model can be written in 
the form 

S + I^2I, (la) 

I^S, (lb) 

where we explicitly have indicated the infection (reaction) rate /3 >0, and the recovery 
(healing) rate ^ > 0. The SIS model is the prototype for a disease spreading model 
where individuals that have recovered from a disease will automatically be susceptible 
to re-infection. On the other hand, if individuals that recover from the disease become 
immune to it, the classic (prototype) epidemic model describing this situation is the 
related susceptible-infected-removed (SIR) model p7) . This latter model, however, 
will not be discussed further in this work. 

The first property to note for the dynamics of the SIS model ([T]l is that it conserves 
the total number of individuals, or in the language of physics, the total number of 
"particles". This number we denote by 

N=Ns{t)+Nj{t), (2) 

where Np{t) is the total number of individuals of type p = I,S at time t. Note that 
in writing Eq. Q, we have explicitly indicated a time-dependence for Np{t) since 
the number of individuals of type p may fluctuate in time even if the total number of 
individuals N in the system remains constant at all times. The second feature of the 
SIS dynamics is the active role played by the infected individuals (the "/ particles"); 
if they have disappeared completely at time to, i.e. Ni{to) = 0, there is no way in the 
model to reintroduce the infection and, therefore, the population remains uninfected 
for all later times t > to. 

2. 1 SIS dynamics in homogeneous space 

In order to understand the dynamics of the SIS model, we will start by studying it in 
homogeneous (continuous) space under the assumptions that susceptible and infected 
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individuals are well mixed. Based on reaction ([T]) we consider the model (in discrete 
time) 

Ns{t + 1) = Ns{t) + ^,Nj{t) - f3 ^''^^^^'^*^ , (3a) 
Nj{t+l) = Nj{t)-f,Nj{t) + (3^^^^^, (3b) 

where the second and third term on the right-hand-side of these equations are as- 
sociated with recovery and infection, respectively. For later convenience, we have in 
writing Eq. (|3]l introduced a characteristic population size, A^x , so that that the effective 
infection rate is 

,_ 13 _ [1;, N«N,; "the small village" 
N^ + N I A, N»N^; "the big city" 

From Eq. Q it follows that in the limit N « N,^, the infection rate /?' is independent 
of the total population size N . This situation is typical for a smaller community where 
one, within a rather short period of time, may be "exposed" to almost the entire 
population of the community. For instance, this may happen during a funeral, a music 
gathering or a Christmas party. For this reason, we will below refer to this type of 
interaction being of "the small village" type. In big cities, on the other hand, the 
probability per unit time to meet, and therefore potentially infect, another individual 
is reduced as the population of the city increases. In this situation we have chosen to 
model the effective infection rate with Eq. Q so that /3' ~ N^^ in the limit where 
N »N^. 



2.2 Stationary state 



Under stationary conditions Nj{t) and Ns{t) are independent of time t, which we 
denote by Nj and Ng, respectively. From Eqs. Q it follows that in a stationary state, 
the following relation must be satisfied 



NsNi 



(5) 



with the additional constraint (see Eq. (|2]i) that Ns + Nj = N. Equation (|5]l is trivially 
fulfilled for Nj = (and Ns = N), meaning that the disease has died out and will 
never reappear. However, Eq. (|5]l is also satisfied if Ns = with the requirement 
that Nj = N - Ns > 0. The (non-trivial) condition reads 



Ni = N 



1 



>0, 



(6) 



which holds if 
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Fig. 1 Schematics of the metapopulation model used in this paper. A real network of subpopulations is 
mapped onto a complex network. Individuals are distributed over the nodes of the network where the 
susceptible (light blue men) and infected (dark) individuals interact on the nodes according to the SIS 
dynamics (Eq. j3j), and travel between the subpubulations (nodes) by diffusion. 

This means that an "epidemic phase" exists (Nj > 0) in the long time Hmit when 
individuals are infected faster than they recover . In all other cases, the only acceptable 
solution to Eq. (|5]l is the non-epidemic solution Nj = and Ns = N. Since the boarder 
between these phases occurs when (3/fi= 1 + N^/N, this point in parameter space is 
referred to as the epidemic threshold p7] . 



3 Migration within a metapopulation 

In the previous section, we considered the dynamics of the SIS model under the 
assumption of an isolated, homogeneous well mixed population. It was shown that 
an epidemic phase may exist. However, humans are not uniformly distributed over 
the surface of the earth. Instead, throughout history, we have had the tendency to 
form local communities ("subpopulations"), i.e. areas where the population density 
is higher than in its surroundings. Today the population density of the world is rather 
heterogeneous, and the communities range in size from small villages, consisting of 
a few tens of people, to mega-cities or urban areas where many millions of people 
are living within a rather small area. A group of such spatially separated populations 
that are coupled to each other in one way or the other is called a metapopulation 
(Fig. [TJ. Within one local community of a metapopulation one may argue that the 
SIS interactions ([T]) may take place essentially within homogeneous space, so that 
the epidemic threshold for that community is given by Eq. (|7]). However, for this 
equation to be correct one is required to neglect the exchange of people between the 
various communities of the metapopulation. Today we travel more than ever before; 
we visit neighboring villages and cities, we travel to further away regional centers, 
and do international travels to far away destinations for business or pleasure. Hence, 



the effect of the ever increasing human travel patterns 1 15 ■ 17 1 is that the population 
centers are coupled, something that in some cases (as we will see below) influences 
significantly the epidemic threshold. 

In our model, the nodes of the network represent population centers and the 
links connecting them correspond to the human traveling routes between such cen- 
ters (Fig.[T]l. An integrated world will be assumed; by this we mean that any node of 
the network can be connected to any other node by following the links of the network. 
Technically this is to say that the network is fully connected and therefore consists of 
one giant component p8). We let the network consist of V nodes (or vertices), and 
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let Np{i,t) denote the number of individuals of type p = I,S at population center 
(or node) i = l,2,...,yat time t. Thus, the population of node i at time t becomes 
N{i,t) = Ng{i,t) + Nj{i,t), and by summing over i, the total p-population be- 
comes Np{t) = T,i Np{i,t). Therefore, the total population of the whole network, 
independent of the type of individual, can be written as the time-independent constant 
N = Ns{t) + Ni{t). Furthermore, the topology of the complex network is described 
by the so-called adj acency matrix W 1 28|29 1 . The components of this matrix, Wij > 0, 



describe the weight or capacity of the traveling route from population center j towards 
population center i. To keep matters simple, we will assume that the links are symmet- 
ric: the capacity of the route j ^ iis equal to that ofi^j (i.e. Wij = Wji). Moreover, 
Wij = means that no direct migration (or travel) is possible from j towards i. 

To obtain and use a detailed and accurate description of the traveling patterns of 
man is a complicated and challenging matter |15| , and it is outside the scope of this 
work. Here we will take a simplistic view and use a stochastic migration model where 
no assumptions are made on how people migrate. Individuals "leaving" population 
center i will choose a route from center i to j with a probability that is proportional 
to the capacity Wji. In the next time step, the same procedure is repeated again. This 
simple migration model is essentially what in physics is called a random walk or 
discrete diffusion process on the complex network | |30||34| . 

If the susceptible and infected individuals travel according to two independent 
diffusive processes, it follows from the conservation of "particles" pT[j33| that the 
pure migration processes (/i = /3 = 0) satisfy 

V 

Np{i,t+1) = Np{i,t) + DpY,T,jNp{j,t) - DpNp{t,t). (8a) 

i=i 

Here Dp denotes the "diffusion constant" for p-type individuals, with < Dp < 1, 
or more precisely, the traveling probability. Moreover, in writing Eq. (j8a]) we have 
defined a transfer matrix, T with elements 



T„ = ^, (8b) 



_ _13 
UJj 



where the total outgoing capacity from node j is denoted 

V 



YW^r (8c) 



Under the assumption /i = /? = 0, the number of S- and /-individuals are separately 
conserved at any time, and the population sizes in the stationary states of the diffusive 
migration process are given by 

Np{i)<^uj„ p=SJ. (9) 

This result is readily obtained by substituting the above expression into Eq. ([8| and 
using that the adjacency matrix is symmetric. The implication of equation (|9| is that the 
size of city /village i in the stationary state is proportional to the total migration capacity 
uji of that population center, and it is independent of what the diffusive constant might 
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be. As a result of Eq. (j9]l, a separate subpopulation size distribution of villages and cities 
will not be put explicitly into the model. Instead such information enters implicitly 
via the distribution of weights uji, at least it does so for the pure diffusion process. 
Equation (|9]) appeals to intuition because larger cities will have transportation systems 
of larger capacity than smaller cities or villages, in order to serve the larger population. 
In passing we note that diffusion on complex networks, as described by Eq. ([8]), has 
previously been used to detect network clusters (communities) |29 , 3T[ - p3| , and for 
the cascading failures of complex networks due to overloading ||34j|, to mention a few 
examples. 



4 A reaction-diffusion approach to epidemic spreading in a metapopulation 

We will now combine the two phenomena from the previous two sections; infection 
(reactions) and migration (diffusion). To facilitate future discussion, we introduce the 
average nodal population density of the network, p, defined as the average number of 
individuals per node (or vertex) 

N 

P=y, (10) 

where it is recalled that V is the number of nodes in the network. In terms of population 
density, Eq. Q becomes 

P = Ps{t) + Pi{t), Pp{t) = ^^, p = S,I. (II) 

The infection and migration processes of individuals in a metapopulation can now be 
simulated, given some initial values for Ns{i,0) and Ni{i,0), a realization of the 
network, and values for the model parameters. The equations on which such simula- 
tions are based are obtained by combining Eqs. ([3]l and (j8]l so that the relevant set of 
equations is of reaction-diffusion type. Motivated by the key question in epidemiology 
studies, we focus on the long time behavior of this reaction-diffusion system. There- 
fore, our primary interest will be the average nodal density of infected individuals 
in the stationary state of the system, pj. The corresponding density of susceptibles is 
readily obtained as pg = p- pj and will thus not be explicitly given. We recall from the 
preceding discussion that a non-trivial result pi > requires that initially Nj{i, 0) > 
for at least one node in the network. Therefore, it will from now onward be assumed 
that Ni{i, 0) > for some i's. Moreover, if Dj is zero, there is no way that an initial 
infection can be transmitted in our model to other parts of the metapopulation. Since 
such a "local infection" is not a very interesting and/or realistic situation, we will 
require that Dj > 0. 



4. 1 Numerical results 



Figure[2]shows simulation results (open symbols) for the stationary state of the (scaled) 
nodal density of infected individuals, pi/p, vs. the average nodal population density 
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(a) Scale-free network, 7 = 3.0 
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Fig. 2 The nodal density of infected individuals in the stationary state (p/) obtained by simulations (open 
symbols) or a mean-field approximation (solid lines). The underlying networks were assumed to have a 
scale-free degree distribution: p(fe) ~ k^'^ with 7 = 3.0 and 7 = 2.5. The parameters used were Dj = 1 
and /3//1 = 2, and the various opens symbols labeled in terms of (V, D5) are: Circles O (10^, 1), squares 
□ (10^, 1), diamonds O (10*, 0) and triangles v (10^, 0). We point out that the curves corresponding to 
the diamonds and tiiangles (black) fall onto the same line since Dg = 0. The simulations were based on a 
(sequential) combination of Eqs. j3) and jsj and the analytic results were obtained from Eqs. {12) and (24|. 



p/Px (with (Ox = N^jV) assuming that = 2. The underlying network connecting 
the subpopulations was assumed to be non-weighted and of scale-free type, p{k) ~ 

, characterized by the exponents 7 = 3.0 (Fig. 2(a) 1 and 7 = 2.5 (Fig. 2(b)| . Non- 
weighted simply means that the population density at given node is proportional its 
number of incoming links (see Eq. (j9]l). The simulation results shown in Fig.|2](as 
solid lines) were obtained under the assumption of treating pp{i,t) as a continuous 
variable to save computation time, and the networks were generated by the MoUoy- 
Reed algorithm [35 1. A stochastic simulation approach similar to that being used in 
Ref. [10 1 was also implemented. It was found to give equivalent results except for a 
higher expenditure of computer time. 

It is clearly seen from Fig. [2] that for sufficiently low nodal population density, p, 
the initial infection will eventually die out and leave the whole population uninfected. 
Under such conditions, there is no potential risk of the initial infection developing into 
a large scale epidemic (pandemic). The same simulation also shows that there exists a 
critical population density, p = pc, above which an average nodal density will result in 
a non-zero fraction of the population being infected even in the long time limit. This 
means that for p > pc, we are in an epidemic -phase defined by pj > 0. Such a behavior 
is characteristic for what in physics is known as a phase-transition [36[, for instance, 
observed when water goes from a solid to a liquid phase with increasing temperature. 

In passing, it should be mentioned that in the simulations performed to produce 
Fig. [2] it was assumed that reaction and diffusion were sequential processes, not simul- 
taneous. This is the same assumption made in the original work |[TOj[TT|. However, 
in a subsequent study [ [37) it was pointed out that the continuous-time limit of the 
discrete process used in [fTOpTl (and in this work) is not well defined if the reaction- 
diffusion process were sequential, and the author proposed a formulation for which a 
continuous-time limit could be obtained. However, doing so renders the formulation 
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more complex with rather similar conclusions. So, we have for the sake of clarity 
assumed a sequential two-step reaction-diffusion process even if we are well aware 
of its shortcoming when taking the continuous limit. 



4.2 Analytical treatment of the model 

The critical average nodal density p = pc is of great importance since it separates 
the non-epidemic from the epidemic phase, or in other words, defines the model's 
epidemic threshold. The purpose of this section is to obtain an analytic expression 
for this critical population density in terms of the model parameters. To this end, we 
again follow pO}[TT] and impose a mean-field approximation where all nodes of the 
same degree k are considered equivalent. We introduce the nodal population density 

Pp.k{t) = —Tr : (12) 

Vk 

where Np^k (t) denotes the total number of p = S, /-type individuals at nodes of degree 
k in the network, and T4 represents the number of nodes of degree k. Just like in the 
simulations we assume a two-step reaction-diffusion process with the reaction step 
being the first of the two. Since the reaction step is local to a node, the corresponding 
equations for the nodal densities pp,k{t) are readily obtained from Eq. (|3]l after a 
division by Vk- The diffusion step, however, needs further discussion because the 
formalism used in Sec. [3] involving transfer matrices, is no longer optimal. For the 
migration of individuals of type p into and away from a node of degree fc, we obtain 
the following nodal density 

Pp,k{t + 1) = Pp,k{t) - DpPpAt) + kY,Pik'\k)^Pp,k'{t), (13) 

k' 

that follows from the conservation of nodal density of individuals of type p = S or: 
I. The physical interpretation of Eq. ( [T3| ) is as follows: the first term on the right- 
hand-side is the density that was akeady present at nodes of degree k at time t\ the 
second term is associated with the diffusion of a fraction Dp of this density away from 
nodes of degree fc; and the last term represents the diffusion into nodes of degree k 
from neighboring nodes of degrees k' . In writing Eq. ( [T3] l, we have introduced the 
conditional probability density, p{k'\k), for a node of degree k being connected to a 
node of degree k' . If the network is assumed to be uncorrected, this probability is 
known to be pS) 

p{k'\k)=p{k')^ (14) 

where p{k) is the degree distribution and (fc) = Y,k kp{k) denotes the average degree. 
By substituting Eq. {T4\ into Eq. ([T3]l and performing the sum over fc' in the resulting 
equation, one is led to 

D fc 

Pp,k{t + 1) = Pp,k{t) - Dppp,k{t) + -j^Ppi^)^ (15a) 
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with 

Pp{t) = Zp(k')Pp.k'{t). (15b) 

k' 

To obtain the final set of coupled equations for ps.k{t) and pi^k{t) we proceed by first 
replacing i by t + 1 in Eq. ( |15a[ ) and then replacing the densities Pp,k{t +1) appearing 
on the right-hand side of the resulting equation by the right-hand side of Eq. ([3]l after 
dividing it by Vk- The result is 

dtps,k{t) = - psAt) + (1 - Ds) [psAt) + wiAt) - Prk{t)] 

+ ^[ps{t) + ppi{t)-pr{t)] , (16a) 

dtPiAt) = - PI At) + (1 - Di) [(1 - fi)piAt) + /?A(t)] 

+ ^[(l-M)p/(<) + mi)] , (16b) 
(k) 

where the reaction kernel 

/,N PsAt)PlAt) . 
Px + Pk 



r{t)=Y,p{k)rk{t), (i6d) 



and its average 



were introduced. In writing Eq. ( [161 ) we have introduced the shorthand notation 
dtps,k = PS At + 2) - PsAt) (and similarly for pj fcfl 

In order to obtain the stationary state of Eqs. ( |16a| l and ( |16b[ ), we require that 
Pp,k{t) = Pp,k (with p = S,I) is independent of time so that 

Ps.k = {l-Ds) [ps.k + m,fc - + ^[PS + f^Pi - f^r] , (17a) 

Pi.k = (1 - Dj) [(1 - ^i)pj^k + + ^ [(1 - + ' (17b) 

where 7^ is given by Eq. ( |16c| i, but with PpAt) replaced by Pp^k and F is defined 
as an average over /\ similar to Eq. ( |16d[ ). Following the strategy from Ref. |To][TT], 
we then multiply equation ( |17a| i by p{k) and sum over k to obtain 

PI = -r, (18) 



' The form of these equations reflects the assumption of a sequential process where a reaction step is 
taking place at every node before a diffusion step. They are similar to those presented in the supplementary 
information accompanying jlO|l Ij with the exception that here a more general reaction kernel has been 
used. 
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which enables us to simphfy Eq. {Tf) to 

ps,k = (1 - Ds) [ps,fc + - + IJ^Ps^ (19a) 

Dik 

Pi,k = (1 - Dj) [(1 - ^l)pJ^k + m] + Tf^P/ . (19b) 

{It) 

An obvious solution to Eq. (19\ is the trivial solution ps,fc = p and p/ ^ = 0. In the 
following we will mainly be interested in non-trivial solutions. 

Equations ([T9| are valid for a general form of the reaction kernel /\. and have 
previously been derived in Refs. [10, 11 1 where it was assumed that the reactions on 
all nodes of the network were either of "small village type" (p « p^) or "big city 
type" (p » px) (and referred to as reaction type 1 and 2, respectively). However, the 
results from Refs. [10,,11| do not apply to a mixed reaction type of the form ( |16c[ ) 
that, depending on the nodal population density, interpolates between type 1 and 2 
reactions. Such dependence on the local density introduces heterogeneous dynamics in 
the network. In the following, we investigate the effect on the epidemic threshold from 
this more general reaction type, which we find more adequate for realistic systems. 



4.3 Two special cases: migrating and non-migrating susceptible individuals 

We will now consider two special cases, namely the situations where the susceptible 
individuals can or cannot diffuse while the infected individuals always are assumed 
to diffuse {Di = 1). For these two cases we will for simplicity consider Ds = and 
-Ds = l. 

Case 1: Non-migrating susceptible individuals (Ds = 0). When the susceptible 
individuals cannot migrate within the metapopulation, it follows from Eqs. (fT9a| that 



/3 f, /3 Ps,kPi,k 
Pi,k = -rk = 1 1 — . (20) 

p p Px+ ps,k + Pl,k 

In the last transition of Eq. pO] ) we have taken advantage of the definition ( |16c| i for 
the reaction kernel. By multiplying Eq. (j20]l by (px + ps.k + Pi,k ) / Pi.k and p{k) (the 
degree distribution), summing over k, and using that p = ps + Pi (Eq- (fTT)), leads to 

Formally this expression is only valid when both P5 > and pj > 0. Moreover, since 
the parameters entering this equation are non-negative, it follows that one must also 
require f3/p> 1. A critical nodal density, p = pc, is found by solving Eq. (j2T]i under the 
assumption of pj = 0. This leads us to write the stationary state of the nodal density 
of infected individuals in the form 




pi=V ^ (22a) 
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where the critical nodal density reads 



Ds = 0, 



(22b) 



which is the same result as for an isolated population (see Eq. ([?])). This means that 
when the susceptible individuals are immobile the epidemic threshold is independent 
of the network's topology. 

Case 2: Migrating susceptible individuals (Dg = 1). We will now address the situ- 
ation where both susceptible and infected individuals are allowed to migrate between 
subpopulations. Under the assumption that Ds = Dj = 1, it follows from Eq. ( [T9] l that 



Ps.k 



(k) 



PS] 



Pl,k 



{k) 



PI, 



(23) 



where the nodal densities of both infected and uninfected individuals are proportional 
to the node degree. The same result is also found for pure diffusion processes where 
no infection takes place at all p2][33| . By introducing relations ( |23| ) into Eq. ([T8|, 
solving the resulting equation for pg and using p= ps + pj, we obtain 



Pi = P- 



Px 



1 + 



fc p 
(fc) P. 



Ds = 1. 



(24) 



This is valid when < pi < p and pi = otherwise. Contrary to the Ds = case, 
Eq. p4l ) shows that pi is a complicated function depending both on local interaction 
parameters (/3, p and p^) and global network parameters {p{k)) when susceptible 
individuals are allowed to migrate. By taking the "small village" and "big city" limits 
of Eq. p4l) one obtains 



PI 



[P 



_ Px (fc)^ 



« 1 



km P 

(k) p. 
(*:) Px 



(25) 



where k^ denotes the maximum degree of the underlying network (note that > 1). 
Notice that the quantity ^ is not only large in the "big city" limit where pj p^ » 1, 
but potentially also in the "small village" limit (p/px « 1) if the network is sufficiently 
heterogeneous so fc„i/ {k) is large enough. 

The epidemic threshold, p = p^, is obtained by setting p/ = in Eq. (|24]i. In this 
case, however, we cannot get pc in closed form. Rather, it is defined implicitly through 



MPx 



1 + 



k Pc 
{k} Px 



(26) 



Since the summation cannot be evaluated explicitly for a general degree-distribution 
p{k), one must resort to numerical methods to calculate pc for a given set of parameters. 
Exceptions do however exist. One simple example is a regular grid, /. e. a homogeneous 
network, where every node has (fc) = fc links. In this case Eq. (|26| yields the same 
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Fig. 3 The epidemic threshold pc (scaled by px), vs. the ratio of the infection and recovery rates, 
for various networks and sizes (V). Equations \22h) or \26) were used to calculate pc when Dg = or 
Dg = 1, respectively. The curve labeled by k con'esponds to a regular grid where all nodes have degree k. 
Likewise, the lines labeled by 7 correspond to a scale-free degree distribution p{k) ~ k^'^ . For all curves 
it was assumed that the infected individuals were mobile with Dj = 1. 



critical density as when the susceptible individuals are immobile (Eq. ( |22b| l), for which 
there is no dependence on the size of the network (V) or on the number of Hnks each 
node has (as long as all nodes have the same number of links). 

Equations ( [22] l and (|24]) are analytic mean-field expressions for the nodal density 
of infected individuals in the stationary state. In Fig.|2]these expressions are compared 
to simulations (discussed in Sec. |4. 1[ ) for a few model parameters; They corroborate 
well with each otheij^] In particular we point out that the epidemic threshold, the point 
at which pi goes from zero to being non-zero, is well described. Even if the agreement 
between the simulation and analytic results is satisfactory, it should be stressed that the 
networks used in the simulations were only characterized by their degree distributions 
and were otherwise random. Specifically, no clustering or community structures were 
present |29| and it is left for future research to investigate how well the mean-field 
results apply to such cases. 



4.4 Epidemic threshold dependence on model parameters 

In the previous section we established that the epidemic threshold pc is well described 
by the mean-field results Eqs. ( |22b| i and ( [26]l. Here we will discuss in more detail how 
Pc depends on the model parameters. In FigJslwe present Pc/px as a function of the rate 

- The agreement between the simulations and mean-field results found in this work is better than what 
has been reported previously in Ref. (10. 11 1. This in particular applies to the position of the epidemic 
threshold and the value of p/ for p > pc- We speculate that the poorer agreement found in Ref. |10|1 Ij 
is caused by too high values for the rates /3 and p being used in the simulations. Since these rates are not 
explicitly given in Ref. ^10|1 Ij we base our speculation on being able to reproduce their reported behavior 
by increasing the values of the rates. 
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ratio PI fjL for migrating (Ds = 1) and non-migrating (D5 = 0) susceptibles for a few 
different networks. The general trend is that the threshold drops as the infection rate 
/3 becomes much larger than the recovery rate /i. Moreover, when Dg = 1 (or more 
precisely Dg > 0), the epidemic threshold also drops when the degree distribution 
of the network becomes heterogeneous (fat-tailed), or when the size of the network 
increases. 

The two overlapping black curves in Fig.[3]correspond to either immobile suscep- 
tibles [Ds = 0) on any network type, or mobile susceptibles [Ds = 1) on a regular 
grid where all nodes have fc = 4 neighbours. The epidemic thresholds for these cases 
are given by the simple expression Pc/px = ~ 1)^ (see Eq. ( |22b| l), and have no 
dependence on network size. 

The two lower curves in Fig.|3]show pc for two scale-free networks {p{k) ~ k^^'^) 
of different size. Within a scale-free metapopulation the local population density varies 
quite dramatically which means that the dynamics in some nodes are of the "small 
village" type and some (typically much fewer) are of the "big city" type, and the rest 
somewhere in between. As the network gets bigger, or the "small village" dynamics 
dominate (px » Pk for most nodes), the epidemic threshold drops. This behavior can 
be deduced analytically from Eq. ( |26) (or Eq. ( [ZSj l) to give|^ 

(27) 

Here we see directly that for a very heterogeneous degree distribution or an infinite 
network {V 00), both of which yields (fc^) » (fc)^, the epidemic threshold tends 
to zero when keeping (3/ /i finite. 

From the general result Eq. (|26]l we can also deduce the relationship between pc 
and px ■ Since this equation is a function of the fraction Pc/Px, we conclude that they are 
proportional to each other The proportionality constant is (/?/ p - 1)"^ (see Eq. ( |22b| i) 
for immobile susceptibles (Ds = 0) and for a regular grid (e.g. fc = 4 neighbours). 
For large and heterogenous networks, or when "small city" dynamics dominate, it is 
given by Eq. ( |27| ). 

5 Summary and conclusions 

In this work we have reviewed and applied a generalized reaction-diffusion approach 
to epidemic spreading within a metapopulation with the overall goal to study the 
impact of city-size heterogeneity. As the epidemic model we assumed the well-studied 
susceptible-infected-susceptible (SIS) system. This model exhibits an epidemic and a 
non-epidemic phase, of which the latter has a nonzero number of infected individuals 
in the long time limit. The transition between these states, the so-called epidemic 
threshold, depends non-trivially on the parameters of the model. These parameters are 
the infection and recovery rates, the cross-over population size separating small village 
dynamics from large city dynamics, and the topology of the network distributing the 
subpopulations within the metapopulation. We envision the network to symbolize 

^ In practice we let jj^ « 1, where km, is the maximum number of links, in Eq. (26). 
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how cities and villages are interconnected in a transportation network, and human 
traveling patterns are assumed as simple as possible; People diffuse randomly. In 
principle, the diffusion rates also enter into the epidemic threshold, but here we only 
addressed either immobile and mobile susceptibles (Ds = or Ds = 1). Instead, we 
paid special attention to a previously unexplored case, namely, city-size dependence 
on epidemic spreading. By city size we mean population density. It is an interesting 
aspect because a disease spreads quite differently in a small city (a village) compared 
to in a metropolitan city, and the most realistic situation is that all city sizes coexist 
within a metapopulation. We therefore formally let the local dynamics in each node 
adapt itself to the instantaneous local population density. In a diffusion model, the 
size of the city is proportional to the number of its links (we assume an undirected, 
non-weighted network). This means that the network's structural features enter the 
dynamics of the disease in two ways. First, how people migrate, and second, the 
dynamics in each node. Previous works have assumed that the individuals are infected 
in the same way in all nodes. Our main result is that we obtain (in the mean-field 
limit) and for a general network an implicit expression for the epidemic threshold 
expressed in terms of the model parameters. This expression shows good agreement 
with numerical simulations and its asymptotic behaviour interpolates between known 
limits. 

Epidemic spreading is a topic of great social importance that also has attracted 
the interest of diverse scientific disciplines for decades. The results presented in this 
paper open for a better understanding of the position of the epidemic threshold and its 
complex interplay between the local infection parameters and the global parameters 
characterizing the network topology for a large class of diseases (SIS). Such knowledge 
may in the future prove critical for reducing the outbreak size, slowing down the speed 
of spreading and/or spatially confining infectious diseases. Moreover, if epidemic 
dynamics were better understood, our contemporary societies could earlier, and at 
better levels of accuracy, predict if a local outbreak has a risk of becoming a global 
pandemic. If an accurate early warning could be given, appropriate measures could 
be taken, for example by adjusting the production of vaccines, with the potential of 
saving numerous lives. 
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